State-of-the-art techniques for calculating spectral functions in 
models for correlated materials 


K. Hallberg, D. J. Garcia, Pablo S. Cornaglia, Jorge I. Facio and Y. Nunez-Fernandez 
Centro Atomico Bariloche and Instituto Balseiro, CNEA and CONICET, 84OO Bariloche, Argentina 


in 



c/:! 

(N 

(N 



u 

c/5 


PACS 71.27.+a - Strongly correlated electron systems; heavy fermions 

PACS 75.40.Mg - Numerical simulation studies 

PACS 71.15. -m - Methods of electronic structure calculations 

Abstract - The dynamical mean field theory (DMFT) has become a standard technique for the 
study of strongly correlated models and materials overcoming some of the limitations of density 
functional approaches based on local approximations. An important step in this method involves 
the calculation of response functions of a multiorbital impurity problem which is related to the 
original model. Recently there has been considerable progress in the development of techniques 
based on the density matrix renormalization group (DMRG) and related matrix product states 
(MPS) implying a substantial improvement to previous methods. In this article we review some of 
the standard algorithms and compare them to the newly developed techniques, showing examples 
for the particular case of the half-filled two-band Hubbard model. 
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0 Introduction. — The research on materials having 
^ strong electron-electron correlations due to interactions in 

'-local orbitals has attracted a great deal of attention in re- 

I cent years. This is due to their fascinating properties like 
J>.high temperature superconductivity, colossal magnetore- 
(N sistance or heavy fermion behavior, and their sensitivity 
00 to external fields which makes them attractive in view of 
applications. In these materials, strong electron correla- 
^^tions play a central role and represent a major challenge 
for the understanding and control of the different phe- 


Os 


allowed for band structure calculations of a large variety 
of correlated materials (for reviews see Refs. , where 

the DMFT accounts more reliably for the local correlations 

. 


nomena. In spite of the important success of the methods 


based on Density Functional Theory (DFT) in the elec- 
tronic structure calculations of weakly correlated materi- 
als, major difficulties are found when dealing with / or d 
• ^ electron systems where the screened local interaction en- 
ergy is of the order of the conduction electron bandwidth. 

^The DFT-based local density approximation (LDA) 
and its generalizations are unable to describe accurately 
the strong electron correlations. 

To overcome these difficulties, the Dynamical Mean- 
Field Theory (DMFT) and its cluster versions were devel¬ 
oped [sJI^ , which allow to extend these methods and treat 
the dynamical electron correlations in a reliable way. The 
DMFT has become one of the basic methods to calculate 
realistic electronic band structure in strongly correlated 
systems. The combination of the DMFT with LDA had 


A recent alternative proposal, the Density Matrix Em¬ 
bedding Theory, DMET, was developed, which relies on 
the embedding of the wave functions of a local cluster 
fragment (instead of the local Green functions) in a self- 
consistent finite environment 13 m and which seems to 
be a good alternative to the DMFT. 

A key point of the DMFT method is the solution of an 
associated quantum impurity problem where the fermionic 
environment of the impurity has to be determined self- 
consistently until convergence of the local Green function 
and the local self-energy is reached. This approach is exact 
for the infinitely coordinated system (infinite dimensions), 
the non-interacting model and in the atomic limit. There¬ 
fore the possibility to obtain reliable DMFT solutions of 
lattice Hamiltonians relies directly on the ability to solve 
(complex) quantum impurity models. 

During the early stages of the development of the 
DMFT, several quantum impurity solvers were proposed 
and used successfully such as the iterated perturbation 
theory (IPT) [^-18 , exact diagonalization (ED) 
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Hirsch-Fye quantum Monte Carlo (HFQMC) Ezf_ 

non-crossing approximations (NCA) [24|, and the numeri- 
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cal renormalization group (NRG) ^5 - 29 . These methods 
were good enough to allow for the calculation of the metal- 
insulator transitions and other low-lying energy proper¬ 
ties. However, they suffered from important drawbacks 
that had to be overcome if other interesting properties 
were to be calculated, such as systems having a larger 
number of bands, other kind of interactions (such as spin- 
orbit and electron-phonon coupling, etc), and interband 
hybridization. Among the drawbacks, one can mention 
the sign problem and the difficulty in reaching low tem¬ 
peratures in the HFQMC algorithm 30 , the failure of 


31). 


to the multiband Gutzwiller wave-function approximation 
[^ . Recently, these low energy techniques have been in¬ 
terfaced with LDA calculations 65-71 . In the sections 


below, we present the RISE formalism an compare some 
results with those calculated within the DMRG frame¬ 
work. 

A recent important extension of the DMFT equations 
concerns its application to treat problems out of equilib¬ 
rium, extending it to the Keldysh formalism 72 -74 . In 


this context, an interesting improvement 75 performs an 


the NCA in obtaining a reliable solution for the metallic 
state, the limitation to few lattice sites, far from the ther¬ 
modynamic limit of the ED and the reduced high-energy 
resolution of the NRG technique (see improvements to this 
in Ref. 


More recently, several other impurity solvers have been 
developed that overcome (at least partially) many of these 
problems among which we can mention the Density Ma¬ 
trix Renormalization Group (DMRG) [32}|3^ , the contin¬ 
uous time quantum Monte Carlo (CTQMC) [ffHiT] and 
the fluctuation exchange approximation (FLEX) [42] . The 
CTQMC works at lower temperatures than the HFQMC 
but still suffers from sign problems for certain models e.g. 
with interband hybridization and, most importantly, it 
also requires an analytical continuation of the Green func¬ 
tions from the imaginary to the real frequency axis which 
makes it unreliable for some physical quantities involving 
higher energy bands 43 . In addition, FLEX is limited 
to a certain range of interaction strengths 44 . We will 
expand on the DMRG variant below. 

Other methods proposed as impurity solvers include 
the equation of motion technique (used in a bath with 
separate low and high energy degrees of freedom for sin¬ 


exact mapping of the action of the equations onto a single 
impurity Anderson model with time-dependent parame¬ 
ters. 

Eor newcommers to the field, it is recommendable to 
visit the TRIQS and ALPS code libraries containing state- 
of-the-art methods for solving interacting quantum sys¬ 
tems 
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The DMRG as an impurity solver of the DMFT. 

— It has been shown that the Density Matrix Renor¬ 
malization Group (DMRG) [TO -81 can be used very reli¬ 


gle and multiple-orbital Hamiltonians) [45} 47 , the quasi- 
continuous-time solver [48| , an improved IPT approach for 
large interactions [^ (which can be compared with the lo¬ 
cal moment approach, also developed to deal with strong 
interactions and arbitrary fillings |50| ), and a two-mode 
approximation based on Gutzwiller variational approach 


51 . Also recently proposed are methods based on exact 
diagonalization (ED) improved by the use of a restricted 
active basis set for the impurity 52 , by a stochastic distri¬ 


bution approach 53 and by an augmented version which 


involves finite temperature and cluster perturbation 54 . 84 - 86 , the block Lanczos approach 87 


ably to solve the related impurity problem within DMFT 
32 - ^ . By using the DMRG to solve the related impurity 
problem, the density of states is obtained directly on the 
real axis (or with a very small imaginary offset) being this 
its major advantage as compared, for example, to QMG 
techniques. In addition, no a priori approximations are 
made and the method provides equally reliable solutions 
for both gapless and gapfull phases. More significantly, it 
provides accurate estimates for the distribution of spectral 
intensities of high frequency features such as the Hubbard 
bands and their structure, which are of main relevance for 
analysis of x-ray photoemission and optical conductivity 
experiments, among others. 

Other techniques using alternative methods for the cal¬ 
culation of dynamical properties within the DMRG have 
been proposed and, more recently, methods using 

the time evolution DMRG algorithm (time evolving block 
decimation, TEED) for the one- and two-orbital mod¬ 
els where shown in [83| (see below). 

In the context of the DMRG and the related ma¬ 
trix product states (MPS) as impurity solvers within the 
DMFT, recent developments include the Kernel Polyno¬ 
mial Method (Ghebyshev expansion for Green functions) 

a pole decompo¬ 


Other promising methods have been proposed based on 
configuration-interaction approximations to ED and from 


the quantum chemical perspective 55 


sition technique within the correction-vector method for 
the dynamics 88 and the application to non-equilibrium 
DMFT using MPS 89 . In this work the authors ex¬ 


In recent years the so-called slave boson approach 61 


in the mean field approximation 62 - 64 has been general 


ized to preserve the symmetries of the Hamiltonian in the 
multiorbital case. The rotationally invariant slave-boson 


mean-held theory (RISE) [M -59 , provides a real-axis de¬ 
scription of the low energy excitations of the system. It’s 
main advantages are the lack of hnite size effects and the 
speed at which solutions of the quantum impurity prob¬ 
lem can be found. The lattice version of RISE is equivalent 


plore other geometries for the impurity bath showing an 
increased efficiency for the star environment. 

It was recently realized that converging the DMFT loop 
on the the imaginary-frequency axis rather than the real- 
frequency axis reduces computational costs by orders of 
magnitude because the bath can be represented in a con¬ 
trolled way with far fewer bath sites and, crucially, the 
imaginary-time evolution does not create entanglement. 
The imaginary time setup can therefore treat much more 
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sophisticated model Hamiltonians, opening the possibil¬ 
ity of studying more complicated and realistic models and 
performing cluster dynamical mean field calculations for 
multiorbital situations. The price to be paid is a reduced 
resolution on the real-frequency axis 
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Model and implementation. — As an example of an 
implementation of the DMRG impurity solver, in this pa¬ 
per we describe the half-filled two-orbital Hubbard model 
on a square lattice including a Hund’s coupling: 


^ ~ ^ , ( 1 ) 


where i, j are the sites of a square lattice and brackets 
indicate nearest neighbors, m indicates each of the two 
orbitals and a is the spin of the electron, whose creation 
and destruction operators are and c, respectively. 

Defining rii and Si as the on-site charge and spin oper¬ 
ators respectively, a rotationally invariant on-site Hamil¬ 
tonian is: 


^ U 2 J 2 

hi — 2 nj 2 ■ 


( 2 ) 


lattice (SL): 


Gsl(w + if]) = ^ 




i^si{kx, ky) + UJ + ir] - T,{uj + irj) 


(4) 

where ujsi{kx,ky) = t{cos{kx) -f cos{ky)) with t = 1/2 to 
have a band of half-width D = 1. All energies are given 
in units of D. The lattice Green function Gsl defines a 
new non-interacting Green function = Gg^ -l-E which 
in turn defines a new hybridization t^A(w) = ur + ig — 
gQ^(uj + irj) which is the seed to restart the cycle. The 
procedure is repeated until converged lattice or impurity 
Green functions are obtained (typically between 5 to 10 
iterations). 


To implement the algorithm we consider (96 97 a gen¬ 


eral representation of the hybridization function in terms 
of continued fractions that define a parametrization of 
A(w -f irj) in terms of a set of real and positive coeffi¬ 
cients. As we are here dealing with two levels, each one 
will be hybridized to its own bath. The hybrization can be 
written simply as a continued fraction (“chain geometry” 

M) 






U) Cllm 


6?. 


(5) 




For ^ = 2U and ferromagnetic J (J< 0) the ground state 
of hi is a triplet and the total Hamiltonian reads: 


H 








^l(T^2cr' 


J Sil ■ Si2 t 


<.ij>m(T 




(3) 


where is the spin operator of orbital m at site i. Ap¬ 
plying DMFT to this model leads to a mapping of the 
original lattice model onto an associated quantum impu¬ 
rity problem in a self-consistent bath. In the particular 
case of the two-orbital Hubbard model, the associated 
impurity problem is the single impurity Anderson model 
(SIAM) with two levels, where the hybridization function 
A(a;), which in the usual SIAM is a flat density of states 
of the conduction electrons, is now to be determined self- 
consistently. 

The DMFT iterations consist of the following: Starting 
from a guessed hybridization A(u}+ig) for the impurity, its 
Green function G{u!+iri) is obtained using some numerical 
method. At this point we introduce the DMRG 32 to cal¬ 
culate the dynamics using Lanczos (921 or the correction 


vectors for an array of discretized energies oj [Mpd . From 
this we can compute the self energy T,{u}+ir]) = G~^ 
where go is the non-interacting Green function correspond¬ 
ing to A(a; -I- ig). The self energy allows us to compute 
the Green function on a lattice, in this case on a square 


The ajm and bjm coefficients represent the on-site en¬ 
ergy and nearest-neighbor hopping for the sites of a non¬ 
interacting chain (not related to the sites of the original 
lattice, Eqs. (1-3)), with Hamiltonian: 

/Nc-l 

HsIAM = Hat -f 'y ) I y ) + 

m, cT \ J = 1 

Nc \ 

+ y ) O'jmCjjnaG'ma + h.C. J (6) 

i = l / 

Here Hat = ~ where N and S are 

the total charge and spin operators of the impurity, dlaa- 
creates an electron with spin proyection cr on orbital m 
of the impurity and creates an electron at site i of 
the non-interacting chain of Nq sites. This cut-off on the 
length of the chains corresponds to a cut-off in the rep¬ 
resentation of the continued fraction for the impurity and 
leads to finite-size effects in the spectra. 

In Fig. [T] we show the DMFT-I-DMRG results for the 
densities of states (DOS), 7rA(w -|- ig) = — ImG(a; -I- ig) 
for several values of the local interaction U and a finite 
ferromagnetic J = —O.I. We have kept around 256 states 
in the DMRG decimation procedure and our results were 
benchmarked with exact diagonalization calculations for 
smaller systems. 

This figure depicts the transition between the insulating 
(large U) and the metallic (small U) regimes, showing the 
lower and upper Hubbard bands in both cases. Since our 
intention in this paper is to briefly review the new methods 
developed to solve the impurity within the DMFT, we will 
not analyze the results in detail, but only mention main 
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Fig. 1: (Color online) Converged densities of states for the 
half-filled two-band Hubbard model on the square lattice 
(1) showing the transition from the metallic to the insu¬ 
lating phases. Here L = 2Nc + 2 = 42 sites, J = —0.1 
and T] = 0.2 (which leads to a slightly enhanced DOS in 
the gap for the insulating states). 


physical properties and the comparison between different, 
albeit complementary, methods. Here we can clearly see 
the existence of some structure in the upper and lower 
Hubbard bands, mainly, the presence of a marked peak in 
the inner, low-energy, side of the bands. This feature is 
ubiquitous in recent results for this model and has yet to 
be fully understood. We can also see an incipient peak in 
the middle of each band in the insulating regime, which 
seems to be reminiscent of the van-Hove peak present in 
the non-interacting half-filled square lattice. In a separate 
and more detailed paper we will show results for the metal- 
insulator transition as a function of J 
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To illustrate other reliable techniques mentioned in the 
text, in Fig. [^we show results using three of the numerical 
methods mentioned above (Chebyshev expansion of MPS, 
NRG and analytically continued QMC) for the calcula¬ 
tion of the spectral density for a related two-orbital Hub¬ 
bard model on the Bethe lattice (see and references 
therein). The observed disagreement at high energies is 
due to different broadening convolutions. 

To complete the comparison between DMFT impurity- 
solvers based on DMRG-related methods, such as MPS, 
we show in Fig. the spectral functions calculated us¬ 
ing other recently developed techniques (for the specific 
model and parameters used, please refer to [^). In this 
case, results based on Fourier transformed real-time dy¬ 
namics with matrix product states (TEBD) show a much 
richer structure than previously used methods based on 
QMG, thus providing a much more reliable tool to study 
electronic structure in correlated materials. 

Rotationally Invariant Slave-Boson Mean Field 
Theory . — In the so-called slave boson approach 



Fig. 2: (Color online) DMFT results for the spectral density 
of a half-filled two-band Hubbard model in the metallic regime 
on the Bethe lattice (extracted from Ref. with permission) 
using three different quantum impurity solvers: Chebyshev- 
expanded matrix product states, numerical renormalization 
group, and continuous-time quantum Monte Carlo. 


[^[^, the Hilbert space of the on-site Hamiltonian hi 
(see Eq. is described using an enlarged space which 
includes fermion (quasiparticle) operators with identical 
structure as the original (physical) operators, and a set of 
auxiliary bosonic degrees of freedom. To recover the orig¬ 
inal Hilbert space, a set of constraints is imposed on the 
fermionic and bosonic degrees of freedom. The advantages 
of this approach become clear in the mean field approxi¬ 
mation, where the bosonic degrees of freedom are replaced 
by their mean values. Since the local interaction terms in 
the Hamiltonian can be represented by a quadratic form 
in the bosons, the mean field approximation reduces the 
fermionic action to a quadratic form in the quasiparticle 
operators, together with a set of parameters (the mean 
value of the bosonic operators) that can be calculated in 
a variational way. The result is a set of non-interacting 
fermions having a reduced intersite hopping matrix ele¬ 
ment due to the interactions. The usual interpretation is 
that the mean field approach produces a description of the 
low energy quasiparticles which have a renormalized mass. 


It was shown by Li et al. 64 studying the single or¬ 


bital Hubbard model that, in order to correctly describe 
spin fluctuations, it is necessary to properly symmetrize 
the bosonic Hilbert space which in turn makes the the¬ 
ory spin-rotational invariant. Recently, these ideas have 
been generalized in a consistent and general way in Ref. 
56 where it was also shown that other advantages, such 


as the ability to study general forms of the interacting 
Hamiltonian, or to describe systems in which the quasi¬ 
particle weight is non diagonal, are in close relation with 
having a rotationally invariant formalism. 


In this approach the auxiliary boson operators have two 
indices {(j^AB} associated to eigenstates |H) and \B) of hi 
having the same charge parity. To obtain a one-to-one 
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Fig. 3: (Color online) Comparison of DMFT spectral fnnc- 
tions of a half-filled two-orbital Hubbard model in the metal¬ 
lic regime on the Bethe lattice (extracted from Ref. [83] 
with permission), using other three different quantum impu¬ 
rity solvers: using real-time dynamics together with MPS, 
QMC-|-Maximum Entropy and QMC -I- Fade approximants. 



Fig. 4: (Color online) Quasiparticle weights Z vs U/Uc{J = 0). 
Top: DMRG results, Uc{J = 0)/D = 3.2. Bottom: RISB 
results for D/T = 10000, Uc{J = 0)/Il = 4.8. 


mapping with the original Hilbert space, time-independent 
Lagrange multipliers Aq and A are used to enforce the 
following constraints: 

^ 4^</.cb(S| 6|A) = d, (7) 

A,B,C 

where O = {1, /I//?} and /I creates a quasiparticle in 
orbital a = ma. 

The physical operators are represented by a linear com¬ 
bination of the quasiparticle operators. 




( 8 ) 


where = (dj,..., d]^), and the R matrix is a function 
of the boson fields such that the matrix elements in the 
new representation are identical to the original ones 
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In the quantum impurity problem, the kinetic energy 
is accounted for by the coupling of the atomic states to 
the electron bath which is described by the hybridization 
matrix and can be written as Hy = idli/^d. Using 
Eq. ([^ this can be written in terms of the quasiparticle 
operators f. We get Hy = f. To construct the action 

of the system, we write the interaction terms as a function 
of the bosons, and the fermionic action, which is quadratic 
on the fermion fields, reads: 




(9) 


where T is the temperature, 


G{iuJn) = [iconi - Vf - A]' 


( 10 ) 


and I the identity matrix. We replace the bosons with 
time-independent complex numbers and integrate the 


fermions, 

Z = Tr [e-^^] (11) 

The Free energy = —Tln[Z] in the saddle point approx¬ 
imation, including the constraints, reads 

Q = -TTrln[-G-i(za;„)] - Aq + ^ {^CdXo 

ACD 

+ 5cDEA-{C\i^M\D))(j)AC ( 12 ) 

where Ea is the eigenenergy of state |A). 

The saddle-point equations are obtained performing 
partial derivatives of Q with respect to the bosons and 
the Lagrange multipliers. 

We have solved the model of Eq. [T] using RISB as an 
impurity solver. In Fig. we compare the results for the 
quasiparticle weight, Z, obtained from the DMRG calcu¬ 
lations and the RISB formalism. It is well known that 
slave bosons methods tend to overestimate the metalicity 
of a system. This is reflected in a larger value of Z for a 
given value of the interaction U, and of the critical inter¬ 
action, 14 , at which the metal-insultar transition occurs. 
The figure shows that in the present problem when U is 
renormalized by the critical interaction of the J = 0 calcu¬ 
lation, a good agreement between the methods is achieved 
for non-zero values of J. Interestingly, for J ^ 0, there is a 
jump from finite values of Z to zero at the metal-insulator 
transition. We used the RISB formalism to analyze the 
behavior of the lattice free energy at the transition and 
confirm that the transition is of first order for finite values 
of J while it is of second order for J = 0. 

Conclusions and perspectives. — Recent advances 
in numerical techniques for the resolution of multior¬ 
bital quantum impurity problems have made possible the 
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implementation of the Dynamical Mean Field Theory 
(DMFT) in model Hamiltonians and Density-Functional- 
based methods for the calculation of strongly correlated 
materials in a realistic way. Here we review several tech¬ 
niques which have proven useful to calculate spectral den¬ 
sities and other physical properties of these materials. 

In particular, we focus on the Density Matrix Renormal¬ 
ization Group (DMRG)-based techniques and compare re¬ 
sults with the Rotationally Invariant Slave Boson (RISE) 
spectra. Both methods are complementary: RISE allows 
for approximate calculations in the thermodynamic limit 
of the low energy properties at a relatively low computa¬ 
tional cost. DMRG (or its equivalent MPS implementa¬ 
tion) is a controlled numerical method which allows for the 
calculation of spectral densities directly on the real axis at 
all energy scales. As an example, we solved a two-orbital 
Hubbard model including a Hand’s coupling term J on a 
square lattice. We found a metal-to-insulator transition as 
a function of the local repulsion U which is second order 
for J = 0 and becomes first order for J ^ 0. These new 
proposals for calculating spectral functions on the real axis 
show a richer structure in the correlated bands when com¬ 
pared to more traditional methods for respponse functions 
on the imaginary axis. 

New techniques relying on the optimization of Hilbert 
spaces from the quantum information perspective 100 


101 such as the DMRG and extended methods includ¬ 
ing MPS and tensor networks in general, are being de¬ 
veloped. They will be able to tackle more complex mul¬ 
tiorbital Hamiltonians and cluster DMFT approximations 
to treat spatially extended correlations. They are likely to 
produce more reliable calculations of higher energy spectra 
with important implications on physical magnitudes such 
as optical conductivity, and non-equilibrium properties of 
complex correlated systems. 
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